	SUBROUTINE SPECTRUM1(X,N,M,R,S)
	REAL(4),DIMENSION(N)::X
	REAL(4),DIMENSION(0:M)::R,B,S
	REAL(4):: XBAR
	INTEGER:: TAO,T
	REAL(4),PARAMETER::PI=3.1415926
	XBAR=0
	DO I=1,N
	  XBAR=XBAR+X(I)
	END DO
	XBAR=XBAR/N
	S1=0
	DO I=1,N
	  S1=S1+(X(I)-XBAR)**2
	END DO
	DO TAO=0,M
	  R(TAO)=0
	  DO T=1,N-TAO
	    R(TAO)=R(TAO)+(X(T)-XBAR)*(X(T+TAO)-XBAR)/S1
	  END DO
	  R(TAO)=R(TAO)/(N-TAO)
	END DO
	DO L=0,M
	  IF(L.EQ.0.OR.L.EQ.M)THEN
	    B(L)=1
	  ELSE
	  B(L)=1./2
	  END IF
	  S(L)=R(0)
	  DO TAO=1,M-1
	    S(L)=S(L)+R(TAO)*(1+COS(PI*TAO/M))*COS(L*PI*TAO/M)
	  END DO
	  S(L)=B(L)*S(L)/M
	END DO
	END
 	!PROGRAM MAIN
 	!INTEGER,PARAMETER::N=144,M=18
 	!REAL(4),DIMENSION(N)::X
 	!REAL(4),DIMENSION(0:M)::R,S
 	!OPEN(10,FILE='SST.DAT')
 	!READ(10,'(12F8.2)')X
 	!CLOSE(10)
 	!CALL SPECTRUM1(X,N,M,R,S) 
 	!OPEN(12,FILE=' SPECTRUM1.DAT')
 	!DO I=0,M
 	!  WRITE(12,' (2F8.3)')R(I),S(I)
 	!END DO
 	!END 

	SUBROUTINE SPECTRUM2(N,M,X1,X2, R12,R21,P11,P22,P12,Q12,R212,THITA)
	REAL(4),DIMENSION(N)::X1,X2
	REAL(4),DIMENSION(0:M)::R12,R21,R11,R22,P12,Q12,P22,P11,R212,THITA,B
	REAL(4):: X1BAR,X2BAR
	INTEGER:: TAO,T
	REAL(4),PARAMETER::PI=3.1415926
	X1BAR=0
	X2BAR=0
	DO I=1,N
	  X1BAR=X1BAR+X1(I)
	  X2BAR=X2BAR+X2(I)
	END DO
	X1BAR=X1BAR/N
	X2BAR=X2BAR/N
	S1=0
	S2=0
	DO I=1,N
	  S1=S1+(X1(I)-X1BAR)**2
	  S2=S2+(X2(I)-X2BAR)**2
	END DO
	S1=SQRT(S1/N)
	S2=SQRT(S2/N)
	DO TAO=0,M
	  R11(TAO)=0
	  R22(TAO)=0
	  R12(TAO)=0
	  R21(TAO)=0
	  DO T=1,N-TAO
	    R11(TAO)=R11(TAO)+(X1(T)-X1BAR)/S1*(X1(T+TAO)-X1BAR)/S1
	    R22(TAO)=R22(TAO)+(X2(T)-X2BAR)/S2*(X2(T+TAO)-X2BAR)/S2
	    R12(TAO)=R12(TAO)+(X1(T)-X1BAR)/S1*(X2(T+TAO)-X2BAR)/S2
	    R21(TAO)=R21(TAO)+(X2(T)-X2BAR)/S2*(X1(T+TAO)-X1BAR)/S1
	  END DO
	  R12(TAO)=R12(TAO)/(N-TAO)
	  R21(TAO)=R21(TAO)/(N-TAO)
	END DO
	DO L=0,M
	  IF(L.EQ.0.OR.L.EQ.M)THEN
	    B(L)=0.5
	  ELSE
	    B(L)=1.
	  END IF
	  P11(L)=R11(0)
	  P22(L)=R22(0)
	  P12(L)=R12(0)
	  Q12(L)=0
	  DO TAO=1,M-1
	    P11(L)=P11(L)+R11(TAO)*(1+COS(PI*TAO/M))*COS(L*PI*TAO/M)
	    P22(L)=P22(L)+R22(TAO)*(1+COS(PI*TAO/M))*COS(L*PI*TAO/M)
	P12(L)=P12(L)+0.5*(R12(TAO)+R21(TAO))*(1+COS(PI*TAO/M))*COS(L*PI*TAO/M)
	Q12(L)=Q12(L)+0.5*(1+COS(PI*TAO/M))*SIN(PI*L*TAO/M)*(R12(TAO)-R21(TAO))
	  END DO
	  P11(L)=B(L)*P11(L)/M
	  P22(L)=B(L)*P22(L)/M
	  P12(L)=B(L)*P12(L)/M
	  Q12(L)=B(L)*Q12(L)/M
	END DO
	DO L=1,M-1
	  R212(L)=(P12(L)**2+Q12(L)**2)/(P11(L)*P22(L))
	  THITA(L)=180/PI*ATAN(Q12(L)/P12(L))
	END DO
	END
 	PROGRAM MAIN
 	PARAMETER(N=120,M=12)
 	REAL(4),DIMENSION(N)::X1,X2
 	REAL(4),DIMENSION(0:M)::R12,R21,P11,P22,P12,Q12,R212,THITA
 	OPEN(10,FILE='AA2.DAT')
 	DO I=1,N
 	READ(10,'(2F8.2)')X1(I),X2(I)
 	END DO
 	CLOSE(10)
 	CALL SPECTRUM2(X1,X2,N,M,R12,R21,P11,P22,P12,Q12,R212,THITA)
 	OPEN(12,FILE=' SPECTRUM2.DAT')
 	WRITE(12,'("  P11  P22  P12  Q12 R212     THITA  ")') 
 	DO I=0,M
 	  WRITE(12,' (6F10.5)')P11(I),P22(I),P12(I),Q12(I),R212(I),THITA(I)
 	END DO
 	END
